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ABSTRACT 

At the highest levels of pulsar timing precision achieved to date, experiments are limited 
by noise intrinsic to the pulsar. This stochastic wideband impulse modulated self-noise limits 
pulsar timing precision by randomly biasing the measured times of arrival and thus increasing 
the root mean square (rms) timing residual. We discuss an improved methodology of removing 
this bias in the measured times of arrival by including information about polarized radiation. 
Observations of PSR J0437— 4715 made over a one-week interval at the Parkes Observatory 
are used to demonstrate a nearly 40 per cent improvement in the rms timing residual with this 
extended analysis. In this way, based on the observations over a 64 MHz bandwidth centred 
at 1341 MHz with integrations over 16.78 s we achieve a 476 ns rms timing residual. In the 
absence of systematic error, these results lead to a predicted rms timing residual of 30 ns in one 
hour integrations; however the data are currently limited by variable Faraday rotation in the 
Earth's ionosphere. The improvement demonstrated in this work provides an opportunity to 
increase the sensitivity in various pulsar timing experiments, for example pulsar timing arrays 
that pursue the detection of the stochastic background of gravitational waves. The fractional 
improvement is highly dependent on the properties of the pulse profile and the stochastic 
wideband impulse modulated self-noise of the pulsar in question. 

Key words: pulsars: general - pulsars: individual (PSR J0437— 4715) - methods:statistical 



1 INTRODUCTION 

Radio pulsars are very stable rotators whose radio emission is re- 
ceived on Earth as series of radio pulses. The measurement of the 
average pulse train time of arrival (ToA) is referred to as pulsar tim- 
ing and leads to the majority of derived pulsar astrophysics. At the 
heart of this methodology lies the comparison of the estimated ToA 
to a model of the pulsar's properties, including astrometric, spin- 
down, and binary (both Keplerian and relativistic) parameters. Mil- 
lisecond pulsars (MSPs), with their short periods, naiTow features 
in their phase-resolved light curves (pulse profiles ) and very low 
spin frequency derivatives, are exceptionally stable ( IMatsakis et al.l 
1 1997 ) and provide the best opportunity for high precision tim- 
ing experiments. A number of authors have achieved timing pre- 
cision of a few hundred nanosecon ds over time-scales of five to 
ten years (e.g.. IVerbiest et aLll2009l and references therein). With 
such high precision, a number of experime nts are possible, such 
as testing theories of relat i vistic gravity (e.g.fHulse & Tavlor 197^; 
IWeisberg & Tavlod llOOSl ; iKramer et aLll2006l) ; detection of plan- 
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ets orbiting a pulsar jWolszczan & Frail|[T992t ISailes etal.]|201lh 
and the most prec ise published mass me asurements of planets in 
the Solar System jChampion et al.ll201(]|); determining properties 
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of th e interstellar medium (ISM, e.g. 
'201 ll); constraining the equation of state of dense matter (e.g. 
Demo rest e t al. 2010); or detecting irregulari ties in terrestria l time- 
scales and developing a pulsa r based one ( iPetit&Tavellalll996l ; 
iRodirj|2008l ; lHobbs et al.ll2012h . 

One particular application of the pulsar timing methodology is 
the detection of gravitational waves, i.e., perturbations of the space- 
time metric predicted by general relativity. These waves cany en- 
ergy and so their emission can be inferred by observing the de- 
crease of the semi-major axis of an emitting binary system of 
the order of 1 c m a day, as detec ted with pulsa r timing o f bina- 
ries (e g..lTavlor & Weisberg 19 82LlStairs et al.ll9 98; Kram er et al.l 
I2OO6I ; Ijacobv et al., ,2006 ; ,Bhat et alj|2008l ; iFreire et al...20"l2h . A 
method of direct detection of gravitational waves relies on their 
impact on the relative distance between Earth and the pulsar with 
a quadrupolar spatial si gnature as the gr avitational wave passes 
through space-time (e.g.. |Jenet et al loam . 

A number of collaborations around the world are observing 
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~ 30 millisecond pulsars in the pursuit of detecting gravitational 
waves with so called pulsar ti ming array expe rime nts. These exper- 
iments , based on early ideas o f lSazhinl digTSh and'p oster & Backed 
( Il990l) . exploit the correlation betwe en the timing residuals for 
pairs of pulsars in different directions dHellings & Downslll983h in 
the sky to measure the spatial signature of the gravitational wave 
background and thus distinguish it from the timing noise, inter- 
stellar medium (ISM) effects, clock or Solar system ephemeris 
errors, which have different multipole signatures or are un corre- 
lated b etween different objects. We note that recently [van Stratenl 
( I2OI3 I) pointed out that instrumental distortions due to imperfect 
polarimetric calibration can corrupt the quadrupolar signature of 
the gravitational waves. 

A number of authors have demonstrated that in order to detect 
gravitational waves, a high prec ision of pulsar tim ing observations 
must be achieved. For example, |jenetetal.l (l2006') showed that 20 
pulsars need to be timed at a precision of 100 to 500 nanoseconds 
over a period of 10 years. While the desired precision has been at- 
tained for a small number of objects, the timing precision of the 
MSPs is worse than the theoretically predicted values for objects 
expected to be the best timers. A prime e xample of this is the c losest 
and brightest MSP, PSR J0437-4715 djohnston et J]|l993h . De- 
spite numerous attempts, the timing precision of this pulsar has 
been significantly worse than expected from the pulse profile and 
sig nal-to-noise rati o. Sandhu et al. ( 1997), van Straten et al. ( 2001) 
and I Verbiest et al" lioosD all achieved progressively better timing 
precision with the latest results being 199 ns rms timing residual 
over 10 years with integration lengths between 10 and 120 min- 
utes, still a factor of 4 or 5 worse than theoretical predictions. An 
important factor in the improvement of the timing precision was re- 
lated to improved methods of dealing w ith polarimetric calibration 
( lBrittonll2000l ; lvan Stratenll2004l2006l) . 

The long standing problem of PSR J0437— 4715 under- 
performing in terms of timing precision is now understood 
(IOslowskierZll201ll. paper I hereafter) which demonstrated that 
the theoretical predictions neglect the contribution of the stochastic 
wideband impulse modulated s elf-noise (S WIMS), a lso referred to 
as pu lse (or phase) jitter (e.g.. lCordes & Sh annon 2010HLiu et al.l 
l20I2h . For pulsars whose peak flux approaches the system equiva- 
lent flux density this contribution cannot be neglected and it biases 
the estimates of the times of arrival (ToA) of the pulse train. Pa- 
per I also presented a method based on decomposition of the resid- 
ual profiles onto eigenvectors of the noise covariance matrix to help 
remove the effects of SWIMS in post processing of timing residu- 
als. We applied this method to PSR J0437-47I5 and achieved a 20 
per cent reduction in the rms timing residual. The current work is 
concerned with improvements to this methodology, allowing even 
better recovery of unbiased ToAs by analysing the polarized flux of 
the pulsar. 

This paper is structured as follows: Section |2] introduces the 
concept of stochastic wideband impulse modulated self-noise and 
its impact on attainable timing precision. Section [3] describes the 
data used and details of processing techniques applied. In Section 
|4]we present an extension to the methodology allowing for removal 
of the bias introduced by SWIMS. The results of applying this 
method to archival observations of PSR J0437— 47 15 are presented 
in Section|5] We discuss these results in detail in Section |6]before 
drawing our conclusions in Section|7] 

The work presented in this paper is a logical extension of 
the work presented in paper I, which should be consulted in or- 
der to fully understand the cuiTent work. In the following section 



we briefly summarise the key concepts of paper I for the reader's 
convenience. 



2 PULSAR TIMING AND THE STOCHASTIC 

WIDEBAND IMPULSE MODULATED SELF-NOISE 

The pulsar timing methodology is based on cross-correlation of the 
observe d pulse profile with a template, typically done in the Fourier 
domain ( lTavloj|l992[) . The derived ToAs and the pulsar model are 
used to form timing residuals as implemented by the TEMP02 soft- 
ware package (Hobbs et al. 2006). As outlined in Section[T] in the 
case of PSR J0437— 4715 the rms timing residual is much larger 
than what is theoretically expected based on the considerations of 
the pulse profile effective width and the signal-to-noise ratio (S/N) 
of the data. In paper I we provided evidence that this is due to the 
presence of SWIMS, i.e., the pulsar itself becomes a source of sig- 
nificant noise that is normally not accounted for in the noise balance 
and predictions of timing precision. 

The amplitude of SWIMS is phase dependent and thus the 
variance of noise is also phase dependent, that is the noise is het- 
eroscedastic. In addition, due to the emission properties of the pul- 
sar, SWIMS exhibits temporal and spectral correlations. Both types 
of coiTelation have adverse consequences. The temporal correlation 
implies that the noise in the pulse profile is not only heteroscedas- 
tic but also not independent between different phase bins. While 
this correlation can amplify the bias introduced by SWIMS into 
the ToA estimation, it allows for removal of this very bias in post- 
processing. Due to the wideband nature of the pulsar intrinsic noise 
the signal at different frequencies is not independent. Therefore the 
ToA estimated from separate parts of the observing radio band are 
correlated thus neutralising the benefit of increasing the instrumen- 
tation bandwidth. Using wideband receivers is still beneficial for 
other reasons as it allows to address many issues related to the inter- 
stellar medium such as dispersion measure variations or scattering 
(e.g.. (You et al.ll2007j) . 

The temporal correlations and heteroscedasticity of SWIMS 
can be fully characterised by the covariance matrix of the pulse 
profiles after subtracting the template profile that was used during 
the ToA estimation. This allows us to study the noise properties 
at any integration time longer than the pulse period as the covari- 
ance matrix scales inversely with the integration length and the rel- 
ative contribution of SWIMS to the covariance matrix is integration 
length independent. 

For a more detailed discussion of SWIMS we refer the reader 
to Section 2 of paper I. 



3 OBSERVATIONS AND DATA PROCESSING 

Dual-polarization data were recorded at the Parkes 64 m radio tele- 
scope using the central beam of the 20 cm multibeam receiver 
( Staveley-Smith et al. 1996) and the se cond generation of Caltech 
Parkes Swinburne baseband recorder dBailesI 120031 : iHotarj l2006l 
CPSR2). A total of ~ 36.5 hours were recorded between the 19th 
and 27th of July 2003 using two bit digitisation with a bandwidth 
of 64 MHz centred at 1341 MHz and split over 128 channels, each 
0.5 MHz wide, resolved into 1024 phase bins, and integrated over 
16.78 seconds. CPSR2 adjusts the sampling thresholds twice every 
minute to maintain linearity of the digitiser response. The incident 
radiation was dedispersed coherently to preserve the narrow fea- 
tures of the pulse profile. The data were processed as in lvan Strateij 
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phase 

Figure 1. A high S/N template of PSR J0437— 4715 in the 20cm observing band formed from 41.7 hours of data. The top left short panel shows the position 
angle swing and is situated directly above a panel showing the total intensity pulse profile in solid line and the Hnear polarization in a dotted line. The top right 
panel shows the Stokes Q parameter of the template, while the bottom left and right panels correspond to Stokes U and V, respectively. The panels on the left 
use the ordinate axes on the left edge while the right hand side panels use the far right ordinate axes. All stokes parameters are on the same flux scale. 



and we refer the reader to this work for details. The Stokes 
parameters in the calibra ted data follow the conventions defined by 
Ivan Straten et al.l feOlOh . Here, we present a brief summary of the 
applied processing. 

As explained in paper I and Section ISTI below, it is beneficial 
to have a large number of sub-integrations available for the analy- 
sis. Therefore, we demonstrate our methodology using a different 
dataset than in paper I. We describe here an archival dataset am- 
ple enough for our purposes. T his dataset can be obtained from the 
CSIRO's Data Access PortaQ faobbs et alj|201ll) . 

We generated a high S/N, calibrated full stokes template pro- 
file using an independent dataset, as follows: 48 observing sessions 
of PSR J0437— 4715, each at least four hours in duration, spread 
over 7.2 years were calibrated using the measurement equation 
modelling technique (van Straten 2004). This technique makes use 
of observations of a pulsed noise diode that injects a polarized refer- 
ence signal into the feed horn. Observations of the radio source 3C 
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218 (assumed to have constant flux and negligible circular polar- 
ization) were used to constrain the boost along the Stokes V axis. 
Then, the five best quality of the 8 hour sessions were integrated 
together to form an initial pulse profile template. The remaining 
43 sessions w ere mat ched to it using the matrix template match- 
ing technique jvan Str aten 2006). A final pulse profile template was 
created from the best 99.4 hours of data processed in the same man- 
ner and integrated hierarchically to minimise quantisation errors 
into a profile with S/N equal to 18.6 x lO"'. The full polarization 
information was retained and is shown in Fig. [T] Throughout this 
publication we use the same template in all contexts. 

The generated template is used for polarimetric calibration of 
the primary dataset spanning a week of observations. Calibration 
was performed using the measurement equation template match- 
ing technique ( van Straten 2013). This method yields superior re- 
sults by using PSR J0437— 4715 as a polarized reference source. It 
combines measurement equation modelling with matrix template 
matching to break the degeneracies described in the Appendix of 
lvanStratenlj2004 . 
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Figure 2. Timing residuals estimated for 16.78 s sub-integrations of data taken during one of the days timed against the standard from Fig.[T] The mean ToA 
estimation error in the whole dataset is 255 ns, whereas the weighted rms of the residual ToAs is 774 ns. The fit has ;^-/d.o.f. of 11.8. For clarity, we have 
plotted the residuals as a function of MJD. 



In both the dataset used for creating the template profile and 
the one used for estimating the covariance matrix, 15 per cent of 
the frequency channels on each side of the band were rejected to 
avoid quantisation noise and aliasing CPSR2 is a 2-bit recorder 
which therefo re introduces significant d istortions into the measured 

i iulse profile I Jenet & Andersonlll998h . A 2-bit correction scheme 
van Stratenl2013h has been applied and the data were searched for 
cases of imperfect correction. We achieve this by measuring the 
X^/d.o.f. of the fit of the pulse profile to a high signal to noise 
ratio (S/N) template, where d.o.f. refers to the number of degrees 
of freedom of the fit. The profiles with worst and uncorrectable 2- 
bit distortions have distinctly higher values of which allows for 
their easy identification and rejection. All of the processing steps 
were per formed using the PSRCHIvj^ data proces sing and analy- 
sis suite iHotan et al.ll2004l : lvan Straten et al.ll2012h . 

The timing residuals for a subset of the 7845 profiles based 
on the 16.78 s sub-integrations used in the analysis are presented 
in Fig. [2] as a function of MJD. The ToAs were estimated using 
the Fourier domain with Markov chain Monte Carlo algorithm 
in PSRCHIVE. The weighted rms timing residual is 774 ns and 
the X^/'i-O-i- equals 11.8. The unweighted rms timing residual is 
slightly higher and equals 806 ns. The mean ToA estimation er- 
ror is 255 ns, and the mean S/N is 302. Based on the findings of 
paper I, we attribute the higher-than-expected rms timing resid- 
ual to SWIMS and the induced ToA estimation bias and x^/d.o.f. 
on the invalid assumptions used in the ToA uncertainty estimation. 
The rms timing residual scales as the square root of the inverse of 
the integration time for pulsars without temporal correlations be- 
tween subpulses and either without long-term timing noise or when 
studied o ver short periods of time (see, e.g., iHelfand et al.lll97^ : 
iLiu et 1; Shannon & Cordes 2012). Based on this extrapo- 

lation, which we have shown to be true for PSR J0437— 4715 (see 



^ In this context by aliasing we mean reflection of the radio frequencies 
outside of the band. 
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Fig. 3 in Paper I) the data would yield an rms timing residual of 
53 ns when derived from 60 minute integrations. 



4 METHOD 

The method described here is a direct extension of what was pre- 
sented in Se ction 4 of paper I and continues to develop the ideas 
presented by'Demorest (2007'). To summarise, this method consists 
of fitting for a scale, phase shift and offset between the observed 
total-intensity pulse profiles and a template profile and calculating 
the covariance matrix of the residual profiles. Next, the eigenvec- 
tors of the covariance matrix are calculated and used as a new basis 
in which the residual profiles are presented. The eigenvectors cor- 
responding to the highest eigenvalues describe the dominant pulse 
shape variations. As a final step, a correlation between the new pre- 
sentation of the residual profiles and timing residuals is found using 
multiple regression and enables a post-processing correction of the 
ToA bias. 

The method presented in paper I yielded a 20 per cent reduc- 
tion of the rms timing residual by exploiting this correlation. How- 
ever, it is unable to correct the bias introduced by the component 
of SWIMS that is correlated with any of the three degrees of free- 
dom that are removed from the data during the template-matching 
fit. Unfortunately, the component of SWIMS correlated with one 
of these degrees of freedom, namely the template time derivative, 
i.e., the phase shift, has a high impact on the bias of the ToA. In 
this work, we extended the previously published method by incor- 
porating the polarized flux of the pulsar. As we demonstrate, this 
allows to recover some of the signal lost in the previous attempt at 
the mitigation of additional ToA scatter introduced by SWIMS. 

The simplest way to implement this extended methodology is 
by considering the pulse profiles of all four stokes parameters, each 
resolved in N\)[„ phase bins, as a 4A'bin profile with only one polar- 
ization. The additional, artificially introduced, phase bins are popu- 
lated with the polarized flux as measured in the stokes parameters. 
That is, the first A'^in phase bins represent the information contained 
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Figure 3. An example of a profile used as input for the methodology presented in this work. The profile is formed from all four stokes parameters spread over 
4A'bin phase bins instead of four seperate profiles with Wtin phase bins each. We used the same profile as in Fig.[T] The input profile is resolved into 1024 pulse 
phase bins and after combing all the stokes parameters into one profile we obtain a pulse profile with 4096 phase bins. 



in the total intensity measurement, the bins from A'^in + 1 to 2A'bin 
represent the Stokes Q measurements, bins from 2A'bin + 1 to 3A'i,„ 
represent the Stokes U flux and the last A'bin phase bins are pop- 
ulated with the circular polarization, i.e., the Stokes V parameter. 
An example of such a profile is shown in Fig. |3] where the mean 
off-pulse values were set equal to zero for each stokes parameter 
independently. 

Before forming the 4A't,in profiles, we performed template fit- 
ting and subtracted the best-fit template from each observation, 
which is necessary to estimate the residual profile covariance ma- 
trix C. As we considered all four stokes parameters, we have six 
free parameters of the fit: a scale factor to account for flux variation 
of the pulsar (intrinsic, related to the effects of interstellar medium, 
or instrumental effects); a phase shift between the pulse profile and 
the template; and four offsets between the pulse profile and the tem- 
plate for each of the stokes parameters. The phase shift and scale 
factor are determined based on the total intensity but applied to all 
four stokes parameters. A unique baseline offset is calculated from 
and applied to each of the stokes parameters. The fitting procedure 
results in six out of 4A'|,in eigenvectors sharing the eigenvalue equal 
zero and form what we refer to as the fit-space. 

After computing the post-template-fit residual we form the 
4A'bin profiles and directly follow the method presented in paper I 
in equations 2 through to 1 10. When exploiting the additional infor- 
mation available in all the stokes parameters, the sizes of some of 
the mathematical quantities involved are different than presented in 
Paper I, namely: the covariance matrices C and D are both 4A'bin by 
4Afbin square matrices; there are 4A'bin eigenvalues associated with 
as many eigenvectors, each of length 4A'(,in; the matrix of projec- 
tions of residual profiles onto eigenvectors has now 4A'bin columns; 
the vector of covariances between the arrival time residuals and the 
projection coefficients ~f has 4A'bin elements; and there are 4A'bin 
values of the Pearson's product moment correlation coefficients 
between the timing residuals and the projection coefficients onto 



one of the eigenvectors. After the analysis is finished the 4A'bin pro- 
files and eigenvectors are converted back to four Wbin profiles, one 
for each of the four stokes parameters. 

Similar to the methodology described in paper I, we have 
performed a number of simulations to ensure that the extended 
methodology does not remove arbitrary phase shifts not related to 
pulse profile deviation from the template profile. We refer the read- 
ers to our previous work for a description of how such simulations 
are realised. Here we report that no removal of arbitrary shift oc- 
curs and thus the method is valid and useful for astrophysical ex- 
periments. 

An implementation of this method is publicly available as part 
of PSRCHIVE as an upgraded version of the previously available 
application "PSRPCA". It requires the GNU Scientific Librar>0 to 
work and can be significantly sped up by using the CULA librar}0 
faumphrev et"ai]|20iq) if an NVIDIA® CUDA™ enabled Graph- 
ics Processing Unit is available. 



5 RESULTS 

Applying our method to the observed dataset of 7845 profiles leads 
to the following: 

• the detection of significant pulse shape variations with at least 
21 significant eigenvectors, 

• a reduction in rms timing residual from 776 ns to 473 ns, a 39 
per cent reduction; and a reduction in ;|^^/d.o.f . from 1 1 .8 to 4.4. 

The ;|^^/d.o.f. remains large as the arrival time uncertainties are 
estimated without accounting for the impact of SWIMS. 

The first three, most significant eigenvectors are plotted in in 
Fig-Ellll and|6] Each panel represents a different stokes parameter, 
from total intensity on top, via Q and U to V at the bottom. We 



Equation 9 in paper I is missing a multiplicative factor ofl/N. 
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Figure 7. Distribution of correlation coefficients between the arrival time 
residuals and projection on to eigenvectors for the actual observations, 
shown by open circles. The threshold level of three robust standard devi- 
ation estimate is marked with dashed lines. The correlation coefficients are 
highly scattered, thus requiring more stringent criteria for selection of the 
significant eigenvectors. 



note that the Stokes Q and U parameters in the first eigenvector 
are highly correlated and anti-correlated with the Stokes U and Q, 
respectively, of the template profile shown in Fig.[T] The Stokes U 
parameter in the first eigenvector has an additional significant peak 
near the leading edge of the total intensity profile as compared to 
the Stokes U of the template profile. The majority of the detected 
pulse profile variance is located in the central parts of the profile. 

We plot the first 150 values of £, in Fig. |7] In Paper I, we 
chose the number of significant eigenvectors by calculating a resis- 
tant and robust estimate of the standard deviation of t, and searched 
for three subsequent values of t, above three times the measured 
standard deviation. Applying the same criteria for automatic selec- 
tion of the number of significant eigenvector results in 301 signifi- 
cant eigenvectors. This value prompted a closer investigation of the 
correlation coefficients and revealed that the automatic selection 
process does not work very well in this case. With a larger num- 
ber of eigenvectors available, more stringent criteria are necessary 
or the number needs to be chosen by investigating the E, values 
manually. In this case, we arrive at a reasonable number of signifi- 
cant eigenvectors by increasing the required number of consecutive 
points above a chosen threshold to five. 



6 DISCUSSION 

We first discuss the details of the extended method and stress the 
differences between this new method and the one presented in pa- 
per I. We then discuss the impact of the extended method on the 
timing precision of PSR J0437— 4715 and on the prospects of pre- 
cision timing. 



6.1 Discussion of the method 

The previously published method used pulse profiles resolved into 
A^bin phase bins but exploited only one stokes parameter. In the ex- 
tension presented in this publication we employ all four stokes pa- 
rameters with the same phase resolution as before. This implies that 
in order for the covariance matrix to be fully determined we now 



need a minimum of AN\,\^ observations. Our conclusions from pa- 
per I still hold and we briefly reiterate them here: a) when the num- 
ber of available observations is lower than 4A'bin> the number of 
eigenvectors with non-null eigenvalues will be limited to the num- 
ber of observations and care needs to be taken when choosing the 
number of significant eigenvectors; b) the number of required ob- 
servations can be reduced by either gating the pulse profile or by 
reducing phase resolution; c) reducing the phase resolution of the 
pulse profile may wash away the temporal correlation of the noise 
which will render the presented methodology less efficient; d) in 
addition, reducing the phase resolution can lead to aliasing har- 
monics of spin frequencies of the observations as, by definition, 
we are dealing with high S/N observations; e) high number of ob- 
servations is highly desirable in order to increase the S/N of the 
covariance matrix estimate and to ensure complete measurement 
of the SWIMS representations in the data; f) given a constant to- 
tal integration time of all available observations, a larger number 
of separate sub-integrations will yield a higher S/N estimate of the 
covariance matrix. 

As in paper I, the eigenvectors are sorted by placing the cor- 
responding eigenvalues in decreasing order such that the greatest 
variance in the data is associated with the first eigenvector and the 
eigenvectors corresponding to the fit-space are the last ones as all 
of the variance in the data in these dimensions is removed during 
the template fitting procedure. Despite the resulting monotonicity 
of the eigenvalues, the regression coefficients used in the ToA esti- 
mate correction scheme are not monotonic as the highest variance 
in profile does not directly translate into the highest bias of the ToA 
estimate. For this reason, applying multiple regression is crucial for 
a successful correction scheme and projections onto more than one 
eigenvector need to be used. 

We note that care must be taken when making any physical 
interpretation of the eigenvectors, even though the detected pulse 
shape variations are likely to be pulsar intrinsic (see Paper I). This 
is due to the fact that our method calculates the eigenvectors of the 
covariance matrix assuming Euclidean geometry, i.e., the measured 
eigenvectors are orthogonal only in a classical sense. The stokes 
parameters span the Minkowski space and to allow any physical 
interpretations of the measured eigenvectors, their orthogonality 
must be ensured in the Minkowski space instead (see, e.g., Renardyj 
Il996h . Four dimensional Euclidian orthogonal transformations may 
not preserve the time-like interval of the stokes parameters. A prime 
example of this effect is visible in Fig.|4l in which the polarized flux 
of the first eigenvector exceeds total intensity. As was the case in 
the method based on only the total intensity, the obtained eigenvec- 
tors are affected by the template matching process, this time with 
six degrees of freedom. This means that if one was to perform a 
simulation with a known pulse profile distortion inserted into data, 
the measured eigenvector would not correspond directly to the in- 
serted noise (see Fig. 6 in paper I and the relevant text). 

The method presented here is superior to the technique pre- 
sented in Paper I as the latter was limited primarily by any pulse 
shape variations correlated with the time derivative of the total in- 
tensity profile. In paper I, the template matching procedure removes 
all signal projected on the fit-space. Variations that correlate with 
the time derivative of the template profile can introduce very sig- 
nificant bias into the ToA estimate without adding much power to 
the residual profile, as can be demonstrated by simulations. By in- 
cluding the information about the statistics of polarized noise, as 
presented in this work, some part of the information previously 
completely removed during the template fitting is retained in the 
Stokes Q, U and V parameters as these three parameters are not 
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Figure 4. First, most significant eigenvector, presented in tlie same way as the template profile in Fig.[T] except we do not show the polarization angle for the 
eigenvector. Note that the flux in total intensity is smaller than flux in the polarized signal. 



used to estimate ToA. Thus, shape variations in these parameters 
correlated with phase shift are not removed. 

The covariance matrix constructed in the way described in 
Section |4] is related to the covariance matri x of the stokes param- 
eters at a give n phase bin as des cribed by Ivan Str aten (2009) in 
equation 28. In Ivan StratenI Jioiol) the author generalised this ma- 
trix to include cross covariances between different polarized states. 
We note that the covariance matrix C used here contains more in- 
formation than either of the two above formulations as it contains 
information both about covariances between the stokes parameters 
at any phase bins as well as between different phase bins. 



6.2 Application to PSR J0437 4715 

For comparison, the method presented in paper I yields timing 
residual with an rms of 642 ns, a 17 per cent reduction; and 
;i^^/d.o.f. of 9.23. The improvement based on this method is some- 
what worse than for the dataset used in paper I, most likely due to 
the lower S/N or residual 2-bit distortions of the observations used 
in this publication. By exploiting the information available in the 
polarization of the incident radiation we are able to improve the re- 
sults of the bias removal from ToA estimates by more than a factor 



of two compared with the previously published approach. We note 
that despite the significant improvement in the attained rms timing 
residual, the results may vary greatly for different pulsars depend- 
ing on their pulse profile and the nature of SWIMS. 

As we are employing the information contained in the polar- 
ized incident radiation, it is useful to compare the rms timing resid- 
ual obtained using the methodology presented in this paper with 
other methods that use not only the total intensity but the other 
three stokes parame ters as well. Two such methods are the matrix 
template matching jvan Strateriir2006) and timing of the invariant 
interval pulse profile ( Brittoij|2000l) . The first allows us to achieve 
higher timing precision by exploiting the fact that for some pulsars 
the polarized signal contains more power in the higher harmonics 
of the Fourier transform of the pulse profile and can compensate 
for some of the systematic eiTors introduced during the calibration. 
The latter minimises the eiTors arising from polarimetry related er- 
rors by calculating the invariant interval and is for example typi- 
cally used for the precision timing o f PSR J0437— 4715 wit hin the 
Parkes Pulsar Timing Array project ([Manchester et al.l20I3l) . 

The ToAs derived using these two methods have an rms value 
of 588 and 771 nanoseconds and x^/^-O-i- of 4.9 and 5.2, respec- 
tively. The value is significantly better than the corresponding 
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value for ToAs derived from total intensity as calculation of the in- 
variant interval significantly reduces the S/N of the data for highly 
polarized pulsars, such as PSR J0437— 4715, thus making the bias 
introduced by SWIMS less obvious by increasing the uncertain- 
ties in the estimated To As. Furthermore, systematic errors might be 
present in the invariant interval based on observations with CPSR2 
due to the automatic gain control used by this baseband system. 

It is possible to correlate the ToA residuals derived from ma- 
trix template matching with the residual profiles expressed in the 
basis spanned by the eigenvectors to obtain an rms timing residual 
of 449 ns and ;(;^/d.o.f. of 2.9. The rms timing residual is similar 
to the value presented in Section|5]and 23.5 per cent lower than the 
input value. The ;f^/d.o.f. value is smaller than presented in Sec- 
tion |5] as the ToA uncertainty estimate based on matrix template 
matching is more realistic because it is the formal least squares un- 
certainty, which takes into account the covariances between the fit 
parameters. 

The method presented here is very sensitive to polarimetric 
calibration of the data and any imperfections in calibration will be 
detectable by presence of characteristic signals in the eigenvectors. 
Any effect that affects the polarization of the pulsar should be de- 
tectable using this method. The fact that the first eigenvector for 
Stokes Q and U are correlated with the Stokes U and Q parame- 



ters in the template, respectively, can be explained by Faraday ro- 
tation the Earth's ionosphere. Based on the International Reference 
lonospher^H (IRI) and International Geomagnetic Reference Fielt^l 
(IGRF) models, we derived a strong correlation with correlation 
coefficient of 0.83 between the relative position angle derived from 
5-min ute integrations using matrix template matching (van Strate^ 
l2006l) and the predicted position angle change induced by the iono- 
sphere. These calculation s were performed using software devel- 
oped by lYan et alj ( I20T l|) and are shown in Fig. [8] Note that the 
ionosphere was quite stable during the week of observations, how- 
ever the calibration solutions were derived at different times of the 
day, corresponding to different angles between the line of sight and 
the Earth's magnetic field. The coincidence of a spike in Stokes I, 
U and V and the mixing of Stokes Q and U is most likely due to 
the projection of the 4A'bin -dimensional space on to the subspace 
that has the six-dimension fit space removed. It does not imply that 
these two types of variation are correlated. 

For pulsar timing experiments, it is often desirable to de- 
rive ToAs based on hour-long observations. However, the tem- 



^ http://iri.gsfc.nasa.gov 

^ http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html 
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poral correlation of the mixing of Stokes Q and U implies that 
the covariance matrix may no longer scale as the inverse of inte- 
gration time. This in turn renders the predictor derived from one 
time-scale inapplicable to other time-scales. In principle, without 
red noise present and with ionosphere related effects removed, we 
should achieve rms timing residual of 30 ns with one-hour inte- 
grations with the dataset presented. Note that this is comparable to 
what m odem observing systems with four times larger bandwidth 
achieve ( [Manchester et alj2013h . To render the predictors useful at 
all scales, calibration solutions need to be derived at higher time 
resolution, or ultimately the state of the ionosphere can be moni- 
tored in near real-time usin g the dual-frequenc y Global Positioning 
System satellites (see e.g., lGamer et al.|[2008h . Both these options 
should allow to remove the mixing of Stokes Q and U and enable 
higher precision experiments. 

PSR J0437-4715 is about 100 times brighter than other 
known MSPs and, as noted in Paper I, unless this object is un- 
usual, we can expect that SWIMS will limit the timing precision 
of other MSPs with the next generation of radio telescopes such as 
the Square Kilometre Array or Five hundred metre Aperture Spher- 
ical Telescope. Based on the expected collecting areas and system 
temperatures of these instruments (Nan 2006; Sch ilizzi et al.l2O10h . 
one can predict the timing precision of future timing experiments. 
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By doing so, iLiu et alj lioTlh concluded that the timing precision 
of MSPs will be greatly improved with these instruments; however, 
it will likely be limited by the impact of SWIMS on the ToA esti- 
mation. This warrants the continued effort into development of ap- 
propriate algorithm to deal with SWIMS, such as the one presented 
in this work. 



7 CONCLUSION 

A solution to the long-standing problem of failing to achieve the 
expected timing precision for PSR J0437— 4715 was presented in 
paper I. In the same work, a method for improving the precision 
of the timing of this pulsar was presented. Here, we address one 
limitation of the previous method and subsequently attain a reduc- 
tion of more than 40 per cent of the rms timing residual, more than 
a factor of two better than the method published in paper I. This 
corresponds to reduction of observing time necessary to achieve a 
given precision by almost a factor of three. Thus this method is 
likely to be very important for many pulsars when observed with 
the next generation of radio telescopes. The rms timing residual 
achieved with the method presented here is lower than can be at- 
tained with any other available method, including matrix template 
matching and timing of invariant interval pulse profiles. Currently, 
the applicability of this method, and thus the timing precision of 
PSR 10437—4715, is limited by ignoring the Earth's ionosphere in 
the data analysis, thus providing a natural next step in the pursuit 
of ultimate timing precision that remains to be addressed in future 
work. After correcting the influence of Earth's ionosphere it will 
be possible to achieve an rms timing residual of about 30 ns in 
one hour integration. Improving the timing precision for a number 
of pulsars will impact a number of astrophysical experiments, in 
particular the detection of gravitational waves with pulsar timing 
arrays. 
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